Set Up

knitr::opts_chunk$set(
    echo = TRUE)

library(tidyverse)
library(tidycensus)
library(sf)
library(kableExtra)
library(tigris)
library(viridis)


options(scipen=999)
options(tigris_class = "sf")

Census API Key

# Make sure to replace this with your own API key! This one is Jenna's.
census_api_key("41e1c0d912341017fa6f36a5da061d3b23de335e", overwrite = TRUE)

Load themes and palettes

mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 16,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    strip.text.x = element_text(size = 14))
}

plotTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 16,colour = "black"),
    plot.subtitle = element_text(face="italic"),
    plot.caption = element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line("grey80", size = 0.1),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    strip.background = element_rect(fill = "grey80", color = "white"),
    strip.text = element_text(size=12),
    axis.title = element_text(size=12),
    axis.text = element_text(size=10),
    plot.background = element_blank(),
    legend.background = element_blank(),
    legend.title = element_text(colour = "black", face = "italic"),
    legend.text = element_text(colour = "black", face = "italic"),
    strip.text.x = element_text(size = 14)
  )
}
qBr <- function(df, variable, rnd) {
  if (missing(rnd)) {
    as.character(quantile(round(df[[variable]],0),
                          c(.01,.2,.4,.6,.8), na.rm=T))
  } else if (rnd == FALSE | rnd == F) {
    as.character(formatC(quantile(df[[variable]]), digits = 3),
                 c(.01,.2,.4,.6,.8), na.rm=T)
  }
}

q5 <- function(variable) {as.factor(ntile(variable, 5))}

# Load hexadecimal color palette
# Discrete from YlGnBu colorbrewer palette
# This way also the lightest is yellow, won't get confused with NA
palette5.YlGnBu <- c("#ffffcc","#a1dab4","#41b6c4","#2c7fb8","#253494")

# Here's a palette of two colors for TOD and Non-TOD, similar to what we used in lab.
paletteTOD <- c("#fdb863", "#80cdc1")

Introduction

Data Wrangling

acs_vars_DC <- c("B02001_001E", # Estimate!!Total population by race -- ##let's double check that it's okay to use this as long as we justify it
                "B02001_002E", # People describing themselves as "white alone"
                "B02001_003E", # People describing themselves as "black" or "african-american" alone
                "B15001_050E", # Females with bachelors degrees
                "B15001_009E", # Males with bachelors degrees
                "B19013_001E", # Median HH income
                "B25058_001E", # Median rent
                "B06012_002E", # Total poverty
                "B08301_001E", # People who have means of transportation to work
                "B08301_002E", # Total people who commute by car, truck, or van
                "B08301_010E", # Total people who commute by public transportation"
                "B03002_012E", # Estimate Total Hispanic or Latino by race
                "B19326_001E") # Median income in past 12 months (inflation-adjusted)



# ---- Washington, DC - Census Data - 2009 ----
tracts2009 <- 
  get_acs(geography = "tract", 
          variables = acs_vars_DC, 
          year=2009, 
          state=11, 
          geometry=T, 
          output="wide") %>%
  st_transform('ESRI:102685') %>% # using the state plane of Maryland for a projection
  rename(TotalPop = B02001_001E, 
         Whites = B02001_002E,
         Blacks = B02001_003E,
         FemaleBachelors = B15001_050E, 
         MaleBachelors = B15001_009E,
         MedHHInc = B19013_001E, 
         MedRent = B25058_001E,
         TotalPoverty = B06012_002E,
         TotalCommute = B08301_001E,
         CarCommute = B08301_002E,
         PubCommute = B08301_010E,
         TotalHispanic = B03002_012E,
         MedInc = B19326_001E) %>%
  dplyr::select(-NAME, -starts_with("B0"), -starts_with("B1"), -starts_with("B2")) %>%
  mutate(pctWhite = (ifelse(TotalPop > 0, Whites / TotalPop,0))*100,
         pctBlack = (ifelse(TotalPop > 0, Blacks / TotalPop,0))*100,
         pctHis = (ifelse(TotalPop >0, TotalHispanic/TotalPop,0))*100,
         pctBlackorHis = (ifelse (TotalPop>0, (Blacks+TotalHispanic)/TotalPop,0)) *100,
         pctBachelors = (ifelse(TotalPop > 0, ((FemaleBachelors + MaleBachelors) / TotalPop),0)) *100,
         pctPoverty = (ifelse(TotalPop > 0, TotalPoverty / TotalPop, 0))*100,
         pctCarCommute = (ifelse(TotalCommute > 0, CarCommute / TotalCommute,0))*100,
         pctPubCommute = (ifelse(TotalCommute > 0, PubCommute / TotalCommute,0))*100,
         year = "2009") %>%
  dplyr::select(-Whites, -Blacks, -FemaleBachelors, -MaleBachelors, -TotalPoverty, -CarCommute, -PubCommute, -TotalCommute, -TotalHispanic)

# ---- Washington, DC - Census Data - 2017 ----

tracts2017 <- 
  get_acs(geography = "tract", 
          variables = acs_vars_DC, 
          year=2017, 
          state=11, 
          geometry=T, 
          output="wide") %>%
  st_transform('ESRI:102685') %>%
  rename(TotalPop = B02001_001E, 
         Whites = B02001_002E,
         Blacks = B02001_003E,
         FemaleBachelors = B15001_050E, 
         MaleBachelors = B15001_009E,
         MedHHInc = B19013_001E, 
         MedRent = B25058_001E,
         TotalPoverty = B06012_002E,
         TotalCommute = B08301_001E,
         CarCommute = B08301_002E,
         PubCommute = B08301_010E,
         TotalHispanic = B03002_012E,
         MedInc = B19326_001E) %>%
  dplyr::select(-NAME, -starts_with("B0"), -starts_with("B1"), -starts_with("B2")) %>%
  mutate(pctWhite = (ifelse(TotalPop > 0, Whites / TotalPop,0))*100,
         pctBlack = (ifelse(TotalPop > 0, Blacks / TotalPop,0))*100,
         pctHis = (ifelse(TotalPop >0, TotalHispanic/TotalPop,0))*100,
         pctBlackorHis = (ifelse (TotalPop>0, (Blacks+TotalHispanic)/TotalPop,0)) *100,
         pctBachelors = (ifelse(TotalPop > 0, ((FemaleBachelors + MaleBachelors) / TotalPop),0)) *100,
         pctPoverty = (ifelse(TotalPop > 0, TotalPoverty / TotalPop, 0))*100,
         pctCarCommute = (ifelse(TotalCommute > 0, CarCommute / TotalCommute,0))*100,
         pctPubCommute = (ifelse(TotalCommute > 0, PubCommute / TotalCommute,0))*100,
         year = "2017") %>%
  dplyr::select(-Whites, -Blacks, -FemaleBachelors, -MaleBachelors, -TotalPoverty, -CarCommute, -PubCommute, -TotalCommute, -TotalHispanic)


# ---- Combining 2009 and 2017 data ----
allTracts <- rbind(tracts2009,tracts2017)

#------Dropping tracts

allTracts <- subset(allTracts, GEOID!="11001006202" & 
                      GEOID!="11001000201" & 
                      GEOID!= "11001008905" & 
                      GEOID!= "11001009809" & 
                      GEOID!= "11001008600" & 
                      GEOID!= "11001006804" & 
                      GEOID!= "11001001600" &
                      GEOID!= "11001000901")

# ---- Wrangling Transit Open Data ----

wmataStops <- 
  rbind(
    st_read("https://opendata.arcgis.com/datasets/54018b7f06b943f2af278bbe415df1de_52.geojson") %>% 
      select(NAME, LINE)) %>%
  st_transform(st_crs(tracts2009))
  

wmataStops <- wmataStops %>%
  mutate(LINE = strsplit(as.character(LINE), ",")) %>%
  unnest(LINE) %>%
  mutate(across(where(is.character), str_trim)) %>%
  mutate(LINE = replace(LINE,LINE == "slvr", "silver"))%>%
  mutate(LINE = replace(LINE,LINE == "yllw", "yellow"))%>%
  mutate(LINE = replace(LINE,LINE == "grn", "green"))
wmataStops

Visualize WMATA

# ---- Visualize wmata stops ----
## removing color for now - we just want ot see where the stops are, not different lines at this point
ggplot() + 
  geom_sf(data=st_union(tracts2009)) +
  geom_sf(data=wmataStops, 
          aes(),
          show.legend = "point", size= 2) +
  ## scale_colour_manual(values = c("red" = "red", "orange" = "orange", "yellow" = "yellow", "green" = "green", "blue" = "blue", "silver" = "gray")) +
  labs(title="WMATA Stops", 
       subtitle="Washington, DC", 
       caption="Source: opendata.dc.gov") +
  mapTheme()

# ---- Bringing in wmata Line data ----
wmataLines <- 
  rbind(
    st_read("https://opendata.arcgis.com/datasets/a29b9dbb2f00459db2b0c3c56faca297_106.geojson") %>%
      select(NAME)) %>%
    st_transform(st_crs(tracts2009))
wmataLines

Creating Buffers

# --- Relating WMATA Stops and Tracts ----

# Creating buffers (in feet - note the CRS) around WMATA stops
# Both a buffer for each stop, and a union of the buffers, and then bind these objects together

wmataBuffers <- 
  rbind(
    st_buffer(wmataStops, 2640) %>% # projection is in feet
      mutate(Legend = "Buffer") %>%
      dplyr::select(Legend),
    st_union(st_buffer(wmataStops, 2640)) %>%
      st_sf() %>%
      mutate(Legend = "Unioned Buffer"))


# ---- Small multiple facet_wrap plot showing both buffers ----

ggplot() +
  geom_sf(data=wmataBuffers) +
  geom_sf(data=wmataStops, show.legend = "point") +
  facet_wrap(~Legend) +
  mapTheme()

# ---- Spatial operations ----

# sf object with ONLY the unioned buffer
buffer <- filter(wmataBuffers, Legend=="Unioned Buffer")

# Showing that we examine the three different spatial selection types 
# to further illustrate why the select by centroids method is best
# based on tracts2009

# Clip
clip <- 
  st_intersection(buffer, tracts2009) %>%
  dplyr::select(TotalPop) %>%
  mutate(Selection_Type = "Clip")

# Spatial selection
selection <- 
  tracts2009[buffer,] %>%
  dplyr::select(TotalPop) %>%
  mutate(Selection_Type = "Spatial Selection")

# Centroid-in-polygon join to see which tracts have their centroid in the buffer
selectCentroids <-
  st_centroid(tracts2009)[buffer,] %>%
  st_drop_geometry() %>%
  left_join(dplyr::select(tracts2009, GEOID)) %>%
  st_sf() %>%
  dplyr::select(TotalPop) %>%
  mutate(Selection_Type = "Select by Centroids")


# Visualizing using Total Population as the fill

allTracts.group.TotalPop <- 
  rbind(clip, selection, selectCentroids)
    
SpatialSelectionTypes2009tracts <-
      ggplot(allTracts.group.TotalPop)+
      geom_sf(data = st_union(tracts2009))+
      geom_sf(aes(fill = TotalPop)) +
      labs(title = "Total Population within 1/2 mi of a WMATA station, 2009", subtitle = "Three types of spatial selection") +
      facet_wrap(~Selection_Type) +
      mapTheme() + 
      theme(plot.title = element_text(size=20))
SpatialSelectionTypes2009tracts

# Centroid joins as above, then "disjoin" to get the ones that *don't* join
# Then  add them all together
# Contains a correct 2009-2017 inflation calculation

allTracts.group <- 
  rbind(
    st_centroid(allTracts)[buffer,] %>%
      st_drop_geometry() %>%
      left_join(allTracts) %>%
      st_sf() %>%
      mutate(TOD = "TOD"),
    st_centroid(allTracts)[buffer, op = st_disjoint] %>%
      st_drop_geometry() %>%
      left_join(allTracts) %>%
      st_sf() %>%
      mutate(TOD = "Non-TOD")) %>%
  mutate(MedRent.inf = ifelse(year == "2009", MedRent * 1.14, MedRent)) 

# Visualize Time/Space Groups

TimeSpaceGrps <-
  ggplot(allTracts.group)+
  geom_sf(data = st_union(tracts2009))+
  geom_sf(aes(fill = TOD)) +
  labs(title = "Time/Space Groups") +
  facet_wrap(~year)+
  scale_fill_manual(values = paletteTOD)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))
TimeSpaceGrps

Four small-multiple (2009 & 2017) visualizations comparing four selected Census variables across time and space (TOD vs. non-TOD)

Visualizing Median Rent

## NEED TO DO: add the wmata line to the legend

mapsMedRent <-
  ggplot(allTracts.group)+
  geom_sf(data = st_union(tracts2009))+
  geom_sf(aes(fill = q5(MedRent.inf)), color = NA, alpha = 0.75) +
  geom_sf(data = buffer, fill = "transparent", color = "red", size = 1.25)+
  geom_sf(data = wmataLines, color = "black", size = 1)+
  scale_fill_manual(values = palette5.YlGnBu,
                    labels = qBr(allTracts.group, "MedRent.inf"),
                    name = "Median Rent ($)\n(Quintile Breaks)",
                    na.value="#f5e9dc") +
  labs(title = "Median Rent by Census Tract, 2009-2017", subtitle = "Real Dollars; Red border denotes areas close to WMATA stations") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))
mapsMedRent

Visualizing Percent of Total Population who Identify as “White Alone”

## NEED TO DO: add the wmata line to the legend

mapsPctWhite <-
  ggplot(allTracts.group)+
  geom_sf(data = st_union(tracts2009))+
  geom_sf(aes(fill=q5(pctWhite)), color = NA, alpha = 0.75) +
  geom_sf(data = buffer, fill = "transparent",color = "red", size = 1.25, linetype = "dashed") +
  geom_sf(data = wmataLines, color = "black", size = 1)+
  scale_fill_manual(values = palette5.YlGnBu,
                    labels = qBr(allTracts.group, "pctWhite"),
                     name = "Percent White") +
  labs(title = "Percent of Population that Identifies as White, 2009 & 2017", subtitle = "Red border denotes areas close to WMATA stations") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))
mapsPctWhite

Visualizing [VARIABLE 3]

Visualizing Percent of People Who Commute by Public Transportation

mapPubCommute <-
  ggplot(allTracts.group)+
  geom_sf(data = st_union(tracts2009))+
  geom_sf(aes(fill = q5(pctPubCommute)), color = NA, alpha= 0.75) +
  geom_sf(data = buffer, fill = "transparent",color = "red", size = 1.25)+
  geom_sf(data = wmataLines, color = "black", size = 1)+
  scale_fill_manual(values = palette5.YlGnBu,
                    labels = qBr(allTracts.group, "pctPubCommute"),
                    name = "Percent Public Transport") +
  labs(title = "Percent of Commuters using Public Transportation, 2009 & 2017", subtitle = "Red border denotes areas close to WMATA stations") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))
mapPubCommute

Table: Comparing Selected Census Variables across Time and Space (TOD vs. non-TOD)

allTracts.Summary <- 
  st_drop_geometry(allTracts.group) %>%
  group_by(year, TOD) %>%
  summarize(Rent = mean(MedRent, na.rm = T),
            Population = mean(TotalPop, na.rm = T),
            Percent_White = mean(pctWhite, na.rm = T),
            Percent_Black = mean(pctBlack, na.rm = T),
            Percent_HispanicLatino = mean(pctHis, na.rm=T),
            Percent_BlackorHispanic = mean (pctBlackorHis, na.rom=T),
            Percent_Bach = mean(pctBachelors, na.rm = T),
            Percent_Poverty = mean(pctPoverty, na.rm = T),
            Percent_CarCommute = mean(pctCarCommute, na.rm = T),
            Percent_PubCommute = mean(pctPubCommute, na.rm = T),
            Median_Income = mean(MedInc, na.rm = T))
            
kable(allTracts.Summary) %>%
  kable_styling() %>%
  footnote(general_title = "\n",
           general = "All Tracts Summary")
year TOD Rent Population Percent_White Percent_Black Percent_HispanicLatino Percent_BlackorHispanic Percent_Bach Percent_Poverty Percent_CarCommute Percent_PubCommute Median_Income
2009 Non-TOD 896.6667 3203.167 27.11839 66.23193 5.736111 71.96804 1.595191 18.67759 49.95842 36.01278 34041.67
2009 TOD 968.2895 3117.295 40.62766 46.73710 9.766179 56.50328 3.249539 17.38676 37.16398 38.71905 39507.26
2017 Non-TOD 1273.5922 3651.913 31.01982 59.77887 8.878615 68.65748 1.730428 19.11664 48.66809 34.14159 40847.66
2017 TOD 1537.9859 3905.056 48.60460 38.16006 11.522585 49.68264 2.984523 15.05316 33.22493 39.31436 52227.48

All Tracts Summary
allTracts.Summary %>%
  unite(year.TOD, year, TOD, sep = ": ", remove = T) %>%
  gather(Variable, Value, -year.TOD) %>%
  mutate(Value = round(Value, 2)) %>%
  spread(year.TOD, Value) %>%
  kable() %>%
  kable_styling() %>%
  footnote(general_title = "\n",
           general = "Summary for All Tracts in Washington, DC")
Variable 2009: Non-TOD 2009: TOD 2017: Non-TOD 2017: TOD
Median_Income 34041.67 39507.26 40847.66 52227.48
Percent_Bach 1.60 3.25 1.73 2.98
Percent_Black 66.23 46.74 59.78 38.16
Percent_BlackorHispanic 71.97 56.50 68.66 49.68
Percent_CarCommute 49.96 37.16 48.67 33.22
Percent_HispanicLatino 5.74 9.77 8.88 11.52
Percent_Poverty 18.68 17.39 19.12 15.05
Percent_PubCommute 36.01 38.72 34.14 39.31
Percent_White 27.12 40.63 31.02 48.60
Population 3203.17 3117.29 3651.91 3905.06
Rent 896.67 968.29 1273.59 1537.99

Summary for All Tracts in Washington, DC
allTracts.Summary
## # A tibble: 4 x 13
## # Groups:   year [2]
##   year  TOD    Rent Population Percent_White Percent_Black Percent_Hispani…
##   <chr> <chr> <dbl>      <dbl>         <dbl>         <dbl>            <dbl>
## 1 2009  Non-…  897.      3203.          27.1          66.2             5.74
## 2 2009  TOD    968.      3117.          40.6          46.7             9.77
## 3 2017  Non-… 1274.      3652.          31.0          59.8             8.88
## 4 2017  TOD   1538.      3905.          48.6          38.2            11.5 
## # … with 6 more variables: Percent_BlackorHispanic <dbl>, Percent_Bach <dbl>,
## #   Percent_Poverty <dbl>, Percent_CarCommute <dbl>, Percent_PubCommute <dbl>,
## #   Median_Income <dbl>

Grouped Bar Graphs: Comparing Selected Census Variables across Time and Space (TOD vs. non-TOD)

allTracts.Summary %>%
  gather(Variable, Value, -year, -TOD) %>%
  ggplot(aes(year, Value, fill = TOD)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~Variable, scales = "free", ncol=3) +
  scale_fill_manual(values = paletteTOD) +
  labs(title = "Indicator differences across time and space") +
  plotTheme() + theme(legend.position="bottom")

Graduated Symbol Maps

intro text

# More wrangling to get centroids for all tracts, both years, filtering by TOD
# some more data wrangling
wmataStopscoord <- wmataStops %>% 
                  st_transform(st_crs('ESRI:102685')) %>% 
                  dplyr::mutate(x = sf::st_coordinates(.)[,1], 
                                y = sf::st_coordinates(.)[,2]) %>%
                  distinct(NAME, .keep_all=TRUE) %>%
                  dplyr::select(-LINE)

wmataBuffersStops <- 
  st_transform(wmataStopscoord, ('ESRI:102685'))%>%
    st_buffer(2640) %>% # projection is in feet
      dplyr::select(NAME, geometry)

joinBuffersStopsTracts <-
  st_join(allTracts.group, wmataBuffersStops)%>%
  filter(!is.na(NAME)) 

# grouping by wmata stop - population
totalPopbyStop <- joinBuffersStopsTracts %>%
  group_by(NAME, year, .add=TRUE) %>%
  summarize(sumPop = sum(TotalPop))%>%
  st_drop_geometry() %>%
  right_join(wmataStopscoord, by="NAME")

Graduated symbol map of population within 1/2 mi of each wmata station

# map
## legend work needed - merge the legend items so sumPop circles are
map.TotalPopbyStop <-
  ggplot()+
  geom_sf(data=allTracts.group, color="white", fill="gray", alpha=0.4)+
  geom_point(data=totalPopbyStop, aes(x=x, y=y, size=sumPop), color='blue', alpha=0.5)+
  scale_size_area(max_size = 15) +
  guides(size=guide_legend("Sum of Population")) +
  geom_sf(data=wmataLines, size=1, color="black", alpha=0.8) + 
  aes() +
  labs(title="Sum of Population within 1/2 mi. of WMATA Stops, by Stop", 
       subtitle="Washington, DC", 
       caption="Data: US Census Bureau; opendata.dc.gov") +
  facet_wrap(~year) +
  mapTheme()
map.TotalPopbyStop

Graduated symbol map of median rent within 1/2 mi of each wmata station

# grouping by wmata stop - medrent
avgMedRentbyStop <- joinBuffersStopsTracts %>%
    group_by(NAME, year, .add=TRUE) %>%
    summarize(AvgRent = mean(MedRent.inf, na.rm = TRUE)) %>%
    st_drop_geometry() %>%
    right_join(wmataStopscoord, by="NAME")

# map
map.avgMedRentbyStop <-
  ggplot()+
  geom_sf(data=allTracts.group, color="white", fill="gray", alpha=0.4)+
  geom_point(data=avgMedRentbyStop, aes(x=x, y=y, color=AvgRent, size=5))+
  scale_size_area()+
  scale_color_viridis(option="D", direction= -1) +
  geom_sf(data=wmataLines, size=1, color="black") + 
  aes() +
  labs(title="Avg. Median Rent within 1/2 mi. of WMATA Stops, by Stop", 
       subtitle="Washington, DC", 
       caption="Data: US Census Bureau; opendata.dc.gov") +
  facet_wrap(~year) +
  mapTheme()
map.avgMedRentbyStop

# grouping by wmata stop - population AND medrent
joinPopRentbyStop <- joinBuffersStopsTracts %>%
  group_by(NAME, year, .add=TRUE) %>%
  summarize(sumPop = sum(TotalPop), AvgRent = mean(MedRent.inf, na.rm = TRUE)) %>%
  st_drop_geometry() %>%
  right_join(wmataStopscoord, by="NAME")

# map
## legend work needed - avg rent bar, put lowest on top, but keep it  yellow
map.gradsymbolcombo <-
  ggplot()+
  geom_sf(data=allTracts.group, color="white", fill="gray", alpha=0.4)+
  geom_point(data=joinPopRentbyStop, aes(x=x, y=y, size=sumPop, color=AvgRent))+
    scale_size_area(max_size = 13) +
    guides(size=guide_legend("Sum of Population"))+
    scale_color_viridis(option="D", direction= -1) +
  geom_sf(data=wmataLines, size=1, color="black") + 
  aes() +
  labs(title="Population and Avg. Median Rent within 1/2 mi. of WMATA Stops, by Stop", 
       subtitle="Washington, DC", 
       caption="Data: US Census Bureau; opendata.dc.gov") +
  facet_wrap(~year) +
  mapTheme()
map.gradsymbolcombo

Crime section

Intro text ## Crime data wrangling

#2009 Crime data wrangling
# 2009
DC_2009_Crime <-
  rbind(
    st_read("https://opendata.arcgis.com/datasets/73cd2f2858714cd1a7e2859f8e6e4de4_33.geojson") %>% 
      select(OBJECTID, OFFENSE)) %>%
  st_transform(st_crs(tracts2009))
## Reading layer `Crime_Incidents_in_2009' from data source `https://opendata.arcgis.com/datasets/73cd2f2858714cd1a7e2859f8e6e4de4_33.geojson' using driver `GeoJSON'
## Simple feature collection with 31248 features and 23 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -77.11276 ymin: 38.81348 xmax: -76.91002 ymax: 38.99422
## geographic CRS: WGS 84
# getting coordinates, 2009
DC_2009_Crime <- DC_2009_Crime %>%
  dplyr::mutate(x = sf::st_coordinates(.)[,1],
                y = sf::st_coordinates(.)[,2])
#2017 Crime data wrangling
# 2017
DC_2017_Crime <- 
  rbind(
    st_read("https://opendata.arcgis.com/datasets/6af5cb8dc38e4bcbac8168b27ee104aa_38.geojson") %>% 
      select(OBJECTID, OFFENSE)) %>%
    st_transform(st_crs(tracts2009))
## Reading layer `Crime_Incidents_in_2017' from data source `https://opendata.arcgis.com/datasets/6af5cb8dc38e4bcbac8168b27ee104aa_38.geojson' using driver `GeoJSON'
## Simple feature collection with 33091 features and 23 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -77.11232 ymin: 38.81348 xmax: -76.91002 ymax: 38.99357
## geographic CRS: WGS 84
# getting coordinates, 2017
DC_2017_Crime <- DC_2017_Crime %>%
  dplyr::mutate(x = sf::st_coordinates(.)[,1],
                y = sf::st_coordinates(.)[,2])
#Some code and some analyses, etc

Multiple ring buffer

# MULTIPLE RING BUFFER -----------------
multipleRingBuffer <- function(inputPolygon, maxDistance, interval) 
{distances <- seq(0, maxDistance, interval)
  distancesCounter <- 2
  numberOfRings <- floor(maxDistance / interval)
  numberOfRingsCounter <- 1
  allRings <- data.frame()
  
  while (numberOfRingsCounter <= numberOfRings) 
  {if(distances[distancesCounter] < 0 & distancesCounter == 2)
    { buffer1 <- st_buffer(inputPolygon, distances[distancesCounter])
      buffer1_ <- st_difference(inputPolygon, buffer1)
      thisRing <- st_cast(buffer1_, "POLYGON")
      thisRing <- as.data.frame(thisRing[,ncol(thisRing)])
      thisRing$distance <- distances[distancesCounter]
    }
    else if(distances[distancesCounter] < 0 & distancesCounter > 2) 
    { buffer1 <- st_buffer(inputPolygon, distances[distancesCounter])
      buffer2 <- st_buffer(inputPolygon, distances[distancesCounter-1])
      thisRing <- st_difference(buffer2,buffer1)
      thisRing <- st_cast(thisRing, "POLYGON")
      thisRing <- as.data.frame(thisRing$geometry)
      thisRing$distance <- distances[distancesCounter]
    }
    
    else 
    { buffer1 <- st_buffer(inputPolygon, distances[distancesCounter])
      buffer1_ <- st_buffer(inputPolygon, distances[distancesCounter-1])
      thisRing <- st_difference(buffer1,buffer1_)
      thisRing <- st_cast(thisRing, "POLYGON")
      thisRing <- as.data.frame(thisRing[,ncol(thisRing)])
      thisRing$distance <- distances[distancesCounter]
    }  
    
    allRings <- rbind(allRings, thisRing)
    distancesCounter <- distancesCounter + 1
    numberOfRingsCounter <- numberOfRingsCounter + 1
  }
  
  allRings <- st_as_sf(allRings)
}

allTracts.rings <-
  st_join(st_centroid(dplyr::select(allTracts.group, GEOID, year)),
          multipleRingBuffer(st_union(wmataStops), 47520, 2640)) %>%
  st_drop_geometry() %>%
  left_join(dplyr::select(allTracts.group, GEOID, MedRent.inf, year),
            by=c("GEOID"="GEOID", "year"="year")) %>%
  st_sf() %>%
  mutate(distance = distance / 5280)

allTracts.rings.group <- allTracts.rings %>%
  group_by(distance,year, .add=TRUE) %>%
  summarise(AvgRent = mean(MedRent.inf, na.rm = TRUE))

ggplot(data=allTracts.rings.group,
       aes(x = distance, y = AvgRent, colour = year)) +
  geom_line(size=2) +
  geom_point(size=3)

## Crime Maps

# Mapping DC crime data 2009
ggplot(subset(DC_2009_Crime, OFFENSE =="ROBBERY")) + 
  geom_sf(data = allTracts.group %>%
          filter(year=="2009"),
          aes(fill = q5(MedRent.inf)), color = NA, alpha=0.75) +
  scale_fill_manual(values = palette5.YlGnBu,
                    labels = qBr(allTracts.group, "MedRent.inf"),
                    name = "Median Rent ($) \n(Quintile Breaks)")+
  geom_point(aes(x=x, y=y), size=0.025) +
  geom_sf(data=wmataStops, size=1.75, shape=22,color="black", fill="#999999")+
  geom_sf(data=buffer, fill="transparent", color="red", size=1.5)+
  labs(title="DC Crime: Robbery Incidents in 2009", 
       subtitle="Washington, DC", 
       caption="Data: US Census Bureau; opendata.dc.gov") +
  mapTheme()

# Mapping DC crime data 2017
ggplot(subset(DC_2017_Crime, OFFENSE =="ROBBERY")) + 
  geom_sf(data = allTracts.group %>%
          filter(year=="2017"),
          aes(fill = q5(MedRent.inf)), color = NA, alpha=0.75) +
  scale_fill_manual(values = palette5.YlGnBu,
                    labels = qBr(allTracts.group, "MedRent.inf"),
                    name = "Median Rent ($) \n(Quintile Breaks)")+
  geom_point(aes(x=x, y=y), size=0.025)+
  geom_sf(data=wmataStops, size=1.75, shape=22,color="black", fill="#999999")+
  geom_sf(data=buffer, fill="transparent", color="red", size=1.5)+
  labs(title="DC Crime: Robbery Incidents in 2017", 
       subtitle="Washington, DC", 
       caption="Data: US Census Bureau; opendata.dc.gov") +
  mapTheme()